import numpy as np
import numpy.linalg as npl
import numpy.random as npr
import numpy.polynomial.legendre as npleg
import newton_raphson as nr


def gradientE(X):
    n=X.size
    result=np.zeros([n])
    for i in range(0,n):
        if (X[i] == -1): #standard choosen : the derivative of |x| takes the value 0 when x=0.
            result[i]=1/(X[i]-1)
        elif (X[i] == 1):
            result[i]=1/(X[i]+1)
        else:
            result[i]=1/(X[i]+1) + 1/(X[i]-1)
        result[i]=result[i]/np.log(10) #the log(10) appears with the derivation of the log.
    return result

def jacobianGradientE(X):
    n=X.size
    result=np.zeros([n,n])
    for i in range(0,n):
        for j in range(0,n):
            if (i==j):
                if (X[i]==1):
                    result[i,j]=-1/((X[i]+1)*(X[i]+1))
                elif (X[i]==-1):
                    result[i,j]=-1/((X[i]-1)*(X[i]-1))
                else:
                    result[i,j]=-(1/((X[i]+1)*(X[i]+1)) + 1/((X[i]-1)*(X[i]-1)))
            else:
                result[i,j]=0
            result[i,j]=result[i,j]/np.log(10)
    return result

def equilibrium(X, N, epsilon):
    return nr.newton_raphson_backtrack(gradientE, jacobianGradientE, X, N, epsilon)

def randomVect(n):
    #Generation of random vectors with coordinates in [-1;1]
    result=np.zeros([n])
    for i in range(0,n):
        result[i]=2*npr.rand() - 1
    return result

def singleComparisonLegendre(X, N, epsilon, epsilon2):
    #This function compares the equilibrium position of the system which initial position are given by n-sized vector X, with the roots of the derivative of the Legendre polynomial of degree n+1. It returns true if  the absolute value of the difference is smaller than epsilon2.
    #The question asked is wether the equilibrium solutions ressemble to the roots of the Legendre polynomials. Let n be the number of position studied. The Legendre Polynomial derivative studied must be a polynomial of degree n because it must have n roots, in order for its roots to be compared with the n the equilibrium position. That is why the Legendre Polynomial used in this function is the one of degree n+1.
    n=X.size
    Y=equilibrium(X, N, epsilon)
    c=np.zeros([n+2]) 
    c[n+1]=1  #c is an array filled with the coefficients of the Legendre polynomial of degree n (i.e. the n+1-sized null vector with a 1 at the end)
    deriv_c=npleg.legder(c) #deriv_c is an array filled with the coefficients of the derivative of the Legendre polynomial of degree n 
    return (npl.norm(Y-npleg.legroots(deriv_c)) < epsilon2)

def globalComparisonLegendre(N, n, epsilon, epsilon2, nb_of_test): 
#N is the maximal number of iterations used in newton_raphson function 
# n is the size of the generated random vectors in the tests.
    tests_OK=0
    tests_KO=0
    for i in range(0,nb_of_test):
        if (singleComparisonLegendre(randomVect(n), N, epsilon, epsilon2)):
            tests_OK+=1
        else:
            tests_KO+=1
    print "Tests OK : "
    print tests_OK
    print "Tests KO : "
    print tests_KO 

#example of test :
# globalComparisonLegendre(100, 3, 0.00001, 0.00001, 100)
